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We present fully ab-initio calculations of van der Waals coefficients for two different situations: 
i) the interaction between hydrogenated silicon clusters; and ii) the interactions between these 
nanostructures and a non metallic surface (a silicon or a silicon carbide surface). The methods used 
are very efficient, and allow the calculation of systems containing hundreds of atoms. The results 
obtained are further analyzed and understood with the help of simple models. These models can 
be of interest for molecular dynamics simulations of silicon nanostructures on surfaces, where they 
can give a very fast yet sufficiently accurate determination of the van der Waals interaction at large 
separations. 



I. INTRODUCTION 



The van der Waals interaction is a common presence 
in the worlds of Physics, Chemistry, and Biology.— Stud- 
ied for more than a century^ it is a dispersive force that, 
being weak at short distances, becomes the dominant at- 
traction between neutral bodies at large separations. It 
results from the non-zero multipole-multipole attraction 
stemming from transient quantum fluctuations. It is the 
interplay between the electrostatic and dispersive inter- 
actions that determines many interesting phenomena in 
Nature, even in the macroscopic world. For example, 
it is van der Waals interactions that are responsible for 
the remarkable ability of geckos to hold to surfaces! 3 - A 
significant example of application of van der Waals in- 
teractions comes from some operating modes of atomic 
force microscopes^ 

However, the realm of van der Waals forces, being 
quantum in nature, is the world of atoms and molecules 
- the nano world. In fact, these forces determine the 
structure of DNA molecules, the folding and dynamics 
of proteins, the adsorption of atoms, molecules or nanos- 
tructures on surfaces, etc. Moreover, the van der Waals 
atom-surface interaction has also been recently studied 
due to their influence on the quantum reflection of ultra- 
cold atoms on surfaces* 5 - In fact, the upsurge of interest on 
Bose-Einstein condensation of ultracold atoms confined 
in magnetic traps^ near a surface should fuel research 
on atom-surface dispersion interactions, since they are 
key factors for the stability of the condensate. They are 
also key ingredients in the building and functioning of 
many of the systems relevant for the emerging fields of 
nanotcchnology and biotechnology. For example, their 
effect was shown to have a profound influence on the os- 
cillatory behavior of microstructures when surfaces are 
in close proximity (100 nm)r£ 



In this context, the purpose of this Article is to present 
fully ab-initio calculations of the van der Waals coeffi- 
cients for the interaction between nanostructures (silicon 
clusters, in particular) and between these nanostructures 
and surfaces. These results are then used to formulate 
simple models that describe with enough precision the 
van der Waals interactions. The models are important, 
as they can be the starting point for, e.g., molecular dy- 
namics simulations of the behavior of silicon nanostruc- 
tures at surfaces. We will be looking at large separa- 
tions, when the overlap between the electronic clouds is 
negligible. At shorter distances, the situation is consider- 
ably more complicated and no satisfying description has 
emerged yet. 

Typically, van der Waals interactions decay with an 
inverse power of the distance between the two bodies un- 
der consideration; The exponent depends on their their 
dimensionality or their metallic character 4 We are inter- 
ested in two specific cases: 

A) The interaction between two finite nanostructures, 
namely atomic clusters of silicon, saturated with hydro- 
gen. We restrict ourselves to the more interesting non- 
retarded regime, i.e. when the time that it takes for the 
photons to travel between the two molecules is negligible. 
In this case, the van der Waals interaction energy has an 
expansion with respect to the inverse of the intermolec- 
ular distance (1/-R) of the form 



oo „ 



(1) 



in terms of the Hamaker constants C n ~£ The first term 
Cq for a pair of molecules A and B, averaged over all pos- 
sible orientations, can be obtained through the relation 
(atomic units will be used hereafter): 
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where o^ x ^ (iu) is the average of the dipole polarizability 
tensor of molecule X, oS x \ evaluated at the imaginary 
frequency iu: 

a ( x )(iu) = iTY[c*W(m)] . (3) 
o 

The higher order terms in the expansion ([T]) can in a sim- 
ilar way be written in terms of higher order polarizability 
tensors. For example, the Cs coefficient will depend on 
the dipole-quadrupole dynamic polarizability (see, e.g., 
Ref. Ilfj ) . In this article we will focus on the leading term 
of the expansion, i.e. the Hamaker constant Cq. 

B) The interaction between a nanostructure and a sur- 
face. Here we will focus on silicon and silicon carbide 
surfaces. In this case, the leading term of the expansion 
of the van der Waals energy as a function of the distance 
Z between the cluster and the surface is proportional to 
Z~ 3 . This fact was first established by Lennard-Jones?ii 
who used a model with a perfectly reflecting metal. The 
theory was later developed by Casimir and Polder^ and 
by LifschitZfi^ For a wide range of particle-substrate dis- 
tances (from approximately 1 nm to even 10 3 nm), the 
interaction energy is given by: 

A ^ = -(i^kF< ^ 

where Zq is a "reference plane" , and C3 is the Lif- 
shitz coefficient. This coefficient can be calculated from 
the dynamical polarizability of the cluster, a(iu), and 
the macroscopic dielectric function of the bulk material, 
6m (iu), both evaluated at imaginary frequencies, iu: 

47r J Q eu{m) + 1 

Note that C3 is expressed only in terms of quantities 
calculated for the bulk crystal. This expression is a gen- 
eral result, also valid for metallic surfaces. The quantity 
that depends on the characteristics of the surface is the 
position of the reference plane Zq. However, for semi- 
conducting surfaces, it can be shownii. that, in absence 
of local field effects, Zq is equal to a/2, where a is the 
interplanar distance. Moreover, it is known that even 
relatively large local field corrections give rise to rather 
small shifts of the reference planed Note, however, the 
position of the reference plane Zq is a more delicate issue 
in the case of a metal, as positioning the reference plane 
at a distance of a/2 from the surface can lead to signifi- 
cant errors in the interaction energy (i.e. about 30% for a 
noble metal surface)^ Further analysis to determine the 
dependence of van der Waals interactions on the surface 
response are under progress. 

It is interesting to remark that it took a long time to 
verify experimentally the predicted Z~ 3 dependence. Al- 
though this term dominates a very wide distance region, 
in the short distance regime it is only a tail of the parti- 
cle potential whose minimum determines the adsorption. 



The first precise measurement of the van der Waals cou- 
plin g b etween an atom and a surface was reported in 
Ref. 15); later on some more experiments have followed. 6 
The experimental difficulties are, in fact, in pair with 
the theoretical ones, since fully ab initio calculations are 
also challenging. Most of the calculations reported in 
the literatur o 17 ^ 18 are generally limited to atoms or very 
small molecules, and the bulk detailed microscopic struc- 
ture is replaced by some model (e.g., the stabilized jel- 
lium model for metals). An interesting approach con- 
sists in modelling the molecular polarizability using the 
form of the long-wavelength density response of a ho- 
mogeneous electron gas. The van der Waals interaction 
is thus efficiently described by the ground-state electron 
densities of the interacting species, obtaining Cq coeffi- 
cients that on average deviate 9% from those obtained 
by TDDFT calculations. 19 All approximations were jus- 
tified by the authors by the difficulty of treating medium 
size molecules within a full ab initio TDDFT approach. 

It is the purpose of our work to demonstrate that cur- 
rent ab initio techniques permit the calculation of the van 
der Waals coefficients for large nanostructures interact- 
ing with realistically described surfaces. The rest of this 
Article is structured as follows: In Section [II] we review 
the methods used to evaluate the dynamical polarizabil- 
ities and the dielectric functions at imaginary frequency; 
in the following section we present the results of our ab 
initio calculations, that are then further analyzed in Sec- 
tion llVI using some simple models. Finally, we draw some 
conclusions in Section [Vj 



II. METHODS 

The main ingredients to evaluate the van der Waals 
coefficients are therefore the electronic polarizability a 
of the cluster and the dielectric constant of the bulk ma- 
terial e, both evaluated at imaginary frequencies. The 
computational methods and the problems involved in the 
calculation of these two quantities are quite different, so 
we will discuss them separately. 



A. Dynamical polarizabilities 

In principle, one can obtain the dynamical polarizabil- 
ities by making use of any quantum-chemistry theory ca- 
pable of handling time-dependent perturbations. As the 
nanostructures we are interested in can be fairly large, 
we choose the time-dependent (TD) extension of density 
functional theory (DFT), since this approach provides an 
excellent compromise between accuracy and feasibility. 
During the past decade, TDDFT— has become one of 
the most important tools to study electronic excitations 
of molecular systems, especially for medium and large 
systems where it is often the only feasible alternative. 

Several different numerical approaches can be found 
in the literature to calculate a at imaginary frequencies 
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within TDDFT. For example, linear response theory can 
be employed to calculate the density-density response 
function \, from which a directly follows: 

aij(u) = J d 3 r J dV n x(r, r', uo) r 3 . (6) 

This approach is quite common and has been used for 
the calculation of the Cg coefficients of molecules and 
clusters. — Alternatively, one can work in real time: By 
propagating the Kohn-Sham equations in real time, it 
is immediate to obtain a(t). A Laplace transformation 
of this quantity yields a(iu) and therefore the Hamaker 
constant Cq . Recently, some of us have shown^ 2 . how this 
procedure can effectively provide Cg coefficients of large 
molecules. A slightly different approach is the polariza- 
tion propagation technique, whose extension to imagi- 
nary frequencies has been used to compute Cq coeffi- 
cients.-— For very large systems, Banerjee and Harbola 
have proposed the use of orbital free TDDFT, providing 
satisfactory results for large sodium clusters^ 

In this article, we use an alternative scheme based on 
the solution of a Sternheimer equation^ It avoids the use 
of empty states, therefore providing a quite good scaling 
(TV 2 ) with the number of atoms. The real-time propa- 
gation technique mentioned above has a similar scaling, 
but we find the prefactor of the Sternheimer approach to 
be smaller. This method has already been used for the 
calculation of many response properties, like atomic vi- 
brations, electron-phonon coupling, magnetic response, 
etc.— In the domain of optical response, it has been 
mainly used for static response, although a few calcu- 
lations at finite (real) frequency have appeared^ 

We have implemented the Sternheimer equation at 
imaginary frequency in the real-space code octopus. 28 
The details are explained in Ref. [29| - although in that 
case the equation is solved for real frequencies. A gen- 
eralization to imaginary frequencies is straightforward. 
The efficiency of the Sternheimer approach is illustrated 
by the size of the clusters studied in this Article: up 
to Sii72H 12 o, i.e., ~300 atoms, computed with relatively 
modest computer systems. Note that, for these systems, 
the time required for the evaluation of the van der Waals 
coefficients is of the same order as the time required for 
the ground-state calculation. 



B. Dielectric constant 

The electronic band structure of bulk semiconductors 
and insulators, which is the starting point to obtain the 
dielectric functions, can nowadays be accurately com- 
puted with ab initio methods^ Much work has been 
done in the past years to determine which approxima- 
tions allow a proper description of electron-electron and 
electron-hole interactions, which is essential to obtain op- 
tical functions (at real frequencies) in agreement with 
experimental data. 31 



The inverse microscopic dielectric function e 1 of a pe- 
riodic system is related to the response function x : 

e- 1 (q, G, G', w) = S G ,& + « (q, G) X (q, G, G', w) , 

(7) 

where q is a vector in the Brillouin zone, G is a reciprocal 
lattice vector, and v is the bare Coulomb interaction. The 
response function \ obeys the matrix equation 

X = Xo + Xo (v + he) X, (8) 

Xo being the independent-particle Kohn-Sham response 
function, and / xc the so-called xc kernel. The macro- 
scopic dielectric function €m can be readily obtained from 
the microscopic e: 

cm (uj) = lim — —, ; r . (9) 

Kl q^o e -i( q ,G = 0,G' = 0,w) KJ 

The simplest approximation that yields the dielectric 
function consists in applying Fermi's golden rule. In 
this approximation, the optical spectrum is calculated 
as a sum of independent transitions between Kohn-Sham 
(KS) or quasiparticle states. This poor man's approach 
is known to exhibit severe shortcomings compared to ex- 
periments.— The next step is the so-called random phase 
approximation (RPA), that includes the effects due to 
the variation of the Hartree potential upon excitation, 
while f xc is set to zero. Unfortunately, the RPA does 
not lead to any significant improvement for most solids, 
especially if there are no particularly pronounced polar- 
izablc inhomogencitics in the charge density. Replacing 
the KS energies with the quasiparticle energies does not 
solve the problem: the peak positions are usually over- 
corrected and the oscillator strength is not modified. 

It is the neglect of variations of the xc potential, which 
include the effect of the electron-hole Coulomb interac- 
tion, that is responsible for an overall disagreement in the 
absorption strength - in particular for the failure to re- 
produce continuum and bound excitons. Unfortunately, 
the adiabatic local density approximation (TDLDA) for 
the xc-kernel in the case of solids is not sufficient to yield 
good dielectric functions. The reason for this failure can 
be traced back to the short-range nature of the TDLDA 
/ xc , while the "exact" / xc is expected to be long-ranged^! 
decaying in momentum space as 1/q 2 . 

A class of kernels that was shown to be yield good re- 
sults is those derived from the Bethe-Salpeter equation 
(BSE) ) 31 i 33 i 34 used together with the quasiparticle band- 
structure. A parameter-free expression, the "nanoquanta 
kernel" , was obtained in several different ways (for a de- 
tailed discussion we refer the reader to Ref. |34j, and refer- 
ences therein) . Although involving a potentially reduced 
computational effort with respect to the BSE, these cal- 
culations are still significantly more cumbersome than 
those within the RPA or the TDLDA. To keep the com- 
putational cost as low as possible, in many cases it is 
enough to use simplified versions of this kernel. It was 
shown that the nanoquanta kernel has the asymptotic 
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form of a long-range contribution (LRC)»i 



/xc atic (q) = - 



35,36 



(10) 



where a statlc is a material dependent parameter, that can 
be related to the dielectric constant. This long-range con- 
tribution alone is sufficient to simulate the strong contin- 
uum exciton effect in the absorption spectrum and in the 
refraction index of several simple semiconductors, like 
bulk silicon or GaAs, provided that quasiparticle ener- 
gies are used as a starting point. A dynamical extension 
of this LRC model^i of the form 



/ x d c yn (q) 



(11) 



leads to remarkable improvements for optical spectra of 
large gap systems with respect to calculations where the 
kernel is imposed to be static. Moreover, the dynami- 
cal approach was proved to be valid also for energies in 
the range of plasmons and for the determination of di- 
electric constants. Note that the parameters of both the 
static and the dynamical model can be related to phys- 
ical quantities, like the experimental dielectric constant 
and the plasmon frequency. 

In this work, we calculated the dielectric functions at 
imaginary frequency using the computer code DP^ an 
ab initio linear response, plane wave, TDDFT code. De- 
spite the enormous amount of studies concerning the ac- 
curacy of different approximations for the xc kernel for 
solids, it is not a priori clear which approximation is more g 
suitable when one wants to work at imaginary frequen- 
cies. In order to clarify this aspect, we tested several 
approximations for the xc kernel. 



III. CALCULATIONS 

A. van der Waals interactions between silicon 
clusters 

We start our discussion by the calculation of Cq be- 
tween the silicon clusters. The clusters were cut from 
bulk silicon, and then saturated with hydrogens along 
the tetrahedral direction of the surface atoms. The ge- 
ometries were optimized with the computer code siesta^ 
employing norm-conserving pseudopotentials, a double ^ 
with polarization basis set, and the PBE parametriza- 
tion^S for the xc potential. 

From the optimized geometries, we then obtained the 
electric polarizability within TDDFT using the Stern- 
heimer equation, as implemented in the computer code 
octopus^ The electron-ion interaction was described 
through norm-conserving pseudopotentials^ and the lo- 
cal density approximation (LDA}*£ was employed in the 
adiabatic approximation for the xc potential. It is known 
that the LDA provides reliable results for semiconducting 
clusters^; Furthermore, from previous experience with 
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FIG. 1: (Color online) Cg B Hamaker constants as a function 
of the square root of the product of silicon atoms in cluster A 
and B. As C§ B scales basically with the product between the 
number of (silicon) atoms in A and B, we plot the Hamaker 
constants divided by this number. Values of Cg B when A and 
B are the same cluster are plotted as red (dark gray) squares, 
otherwise they are plotted as orange (light gray) dots. 




FIG. 2: (Color online) London effective frequency uii - blue 
(dark gray) crosses, and static polarizabilities per atom - red 
(light gray) squares, for the silicon clusters under study, as a 
function of the number of silicon atoms. 



optical spectra calculations^ we know that these results 
will not change significantly with the use of the more so- 
phisticated GGAs. The equations, in this code, are rep- 
resented in a real-space regular grid, whose spacing is 
chosen to be 0.275 A. The simulation box is constructed 
by joining spheres of radius 4.5 A, centered around each 
atom. The integrals in Eqs. ([2|) and |5j) were performed 
with a Gauss-Legendre quadrature using 6 frequency val- 
ues. With these parameters, we estimate the accuracy of 
our numerical calculations to better than 5%. 

In Fig. Q]we show our results for the Gg Hamaker con- 
stant between silicon clusters. As the value of C6 scales 
with the product of the atoms in the cluster A (-/V<^) and 
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FIG. 3: (Color online) Van der Waals C3 coefficients be- 
tween silicon nanoclusters and a silicon surface. The C3 co- 
efficients were divided by the number of silicon atoms in the 
cluster. The different curves were calculated using different 
approximations for the dielectric constant of the bulk crystal 
at imaginary frequencies (see the text for details). 



in cluster B (N^), we divided the Hamaker constant by 
Nj^N<^ to eliminate this dependence. We show both con- 
stants between two identical clusters (homo-molecular - 
as red (dark gray) squares, also presented in Table |TJ 
and between different clusters (hetero-molecular - or- 
ange (light gray) dots). We see that the largest C6 per 
atom squared comes from the interaction between two 
SiH4 clusters, and then the values decrease rapidly until 
slightly above 220 a.u, where it saturates. The few clus- 
ters that fall far from the line are the most asymmetric, 
for which a description in terms of the average of the 
dipole polarizability tensor is not necessarily as good. 

The polarizability at imaginary frequencies can be 
modelled in the London approximation by introducing 
two adjustable parameters: the static polarizability cv(0) 
and one effective frequency lu±: 



a(iu) 



a(0) 



1 + 



(12) 



If we insert {HJl in © we obtain a simplified expression 
for the homo-molecular Hamaker constant in terms of 
these parameters: 



c 6 = ^V(o). 



(13) 



As we have calculated both Cq and a(0) within TDLDA, 
it is easy to extract u>\ from Eq. (fT3| . The resulting effec- 
tive frequencies are plotted in Fig. [21 together with the 
calculated static polarizabilities per number of Si atoms. 
We can observe that u>\ decreases with the number of Si 
atoms, but the dependence on the size of the cluster is 
rather weak, except for the singular case of the smallest 
aggregates. 



FIG. 4: (Color online) Van der Waals C3 coefficient between 
silicon nanoclusters and a silicon carbide surface. The mean- 
ing of the curves is as in Fig. [3] 
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FIG. 5: (Color online) Dielectric function of Si along the imag- 
inary axis as a function of the frequency in a logarithmic scale. 
Results obtained with different approximations for the xc ker- 
nel are compared to the curve extracted from experimental 
data by using Eq. (TTiT l 



van der Waals interactions between silicon 
clusters and dielectric surfaces 



Next, we consider the case of a silicon cluster in prox- 
imity of a surface of Si or SiC in the zincblende phase. 
In this case we want to calculate C3 coefficients, that are 
determined both by the dynamical polarizability of the 
cluster and the dielectric function of the bulk crystal. 

The ground state calculations for the bulk crystals were 
performed using the plane-wave code ABINIT— with 
norm-conserving Hamann pseudopotentials^ for Si and 
C. We used a cutoff energy for the plane wave basis of 
12.5 Ha for Si and 30 Ha for SiC. The unit cell was relaxed 
within the LDA approximation, yielding lattice param- 
eters with an error smaller than 3%. The Kohn-Sham 
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FIG. 6: (Color online) Dielectric function of SiC along the 
imaginary axis as a function of the frequency in a logarithmic 
scale. We compare results obtained within different approxi- 
mations for the xc kernels. — 

energies and wavefunctions yielded by ground state cal- 
culations were employed to calculate dielectric functions 
at imaginary frequencies using the code DPi^£ For the 
response calculations a shifted k-point grid of 256 points 
was used both for Si and SiC. More detailed information 
on the numerics and convergence issues can be found in 
Ref. M 

As discussed before, previous tests on the effect of dif- 
ferent approximations for the xc potential demonstrate 
that the dynamical polarizability of the hydrogenated Si 
clusters is accurately described within the TDLDA, and 
therefore the C§ coefficients are not going to change by 
more than 5% by using different approximations. We de- 
cide thus to focus on the effect of different models for the 
xc-kernel in the calculation of C3 . 

In Figs. and 2] we can compare C3 coefficients for Si 
clusters on a Si surface and Si clusters on a SiC surface, 
respectively. We present results obtained within the RPA 
(violet upright triangles) , the TDLDA (beige inverted tri- 
angles) , using the static LRC kernel (blue diamonds) and 
the dynamical LRC kernel (green crosses). Note that, in 
the case of the RPA and the TDLDA calculations, we 
used the Kohn-Sham band structure to build xoj while 
the GW quasiparticle states are used when the static or 
dynamical LRC kernels are employed (see Refs. |36| and 
1371 for details). In Si and SiC the GW corrections to the 
bandstructures are essentially equivalent to a rigid shift 
of the conduction states, thus replacing KS energies with 
quasiparticle energies leads to a rigid shift of the absorp- 
tion spectrum towards higher energies. We also plotted 
in Figs. [3] and 2] the values of C3 obtained using simple 
models (red squares) for both the dynamical polarizabil- 
ity of the cluster and the dielectric function of the crystal 
at imaginary frequencies. We will discuss these analytical 
models and the quality of their results in Sect ion HVl The 
peaks of C 3 as a function of the number of Si atoms oc- 
cur for higher polarizable clusters. The oscillations of the 



polarizabilities in turn exactly correlate with the binding 
energy, with largest polarizabilities corresponding to the 
most stable clusters. 

For the interaction of Si clusters on either a Si or a 
SiC surface, all the approximations used for the xc ker- 
nel give curves with very similar trends and a dispersion 
of the values which is smaller than 10%. This finding 
reflects the fact that the dielectric function at imaginary 
frequency is a very smooth and well behaved curve, and 
therefore fairly simple to reproduce; it starts at the value 
of the static dielectric constant at iu = 0, and then de- 
creases monotonically to its asymptotic limit of one (see 
Figs. M and [6]). 

Dielectric constants in the imaginary frequency 
axis can be obtained experimentally by performing a 
Kramers-Kronig transformation of the values obtained 
in the real axis, 

e M (iu) = l + - dcu n 14 

as long as the experimental absorption spectra has been 
measured on a large enough spectral range. We include 
in Fig.[5]the experimental curve for Si?i£ for comparison. 
The curve calculated with the dynamical LRC approxi- 
mation is exactly superposed to the experimental curve. 
In fact, this is the only approximation that yields a good 
dielectric constant, which fixes the interception with the 
y axis, and an overall good shape of the absorption spec- 
trum over a large spectral ranged 

It is interesting to notice that the static LRC results 
are worse than even the RPA curve, being overestimated 
over the whole frequency range. In the RPA, there is a 
"fortuitous" compensation of errors (the error due to the 
too high dielectric constant is balanced by the shift of the 
spectral weight to lower energies due to the DFT-LDA 
underestimation of the absorption edge). The TDLDA 
curve is the one that lays further from the dynamical 
LRC solution at lower frequencies due to the even higher 
dielectric constant, but it greatly improves for u > 0.1, 
thanks to the same compensation of errors already ob- 
served for the RPA calculation. 

The same conclusions can be obtained for SiC (see 
Fig. [6]). In this case, we do not have access to exper- 
imental results, but it is reasonable to expect that the 
dynamical LRC approximation will yield the most accu- 
rate result overall. The static LRC results again shows a 
consistent overestimation of the dielectric function, while 
the the RPA curve is the closest to the dynamical LRC 
result. 

In the light of this analysis, one can interpret the re- 
sults for the C3 coefficients. First of all, the dynamical 
LRC kernel is expected to work very well. Outside the 
limits of validity of the dynamical LRC model (large gap 
insulators, strongly bound excitons) only a calculation 
for the bulk crystal based on the solution of the BSE (or, 
equivalently, based on the fully ab initio Bethe-Salpeter 
derived kernel) can guarantee the quality of the C3 co- 
efficients. This implies necessarily larger computational 
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costs. Perhaps surprisingly, the RPA and TDLDA ap- 
pear to be good approximations to evaluate C3, despite 
their well known deficiencies in the calculation of optical 
absorption spectra. This is not necessarily true for every 
system, but it is probably true provided that the calcu- 
lated dielectric function at zero frequency is larger than 
the experimental one. In this case, in fact, we can expect 
the (at least partial) cancellation of error between the 
too high starting point of the curve e(kt) and its too fast 
decay to one as a consequence of the shift of the spec- 
tral weight to lower energies. One should be very careful 
not to use the static LRC approximation for the kernel, 
despite the fact that it gives an absorption spectrum in 
overall agreement with the experiment. This is due to the 
fact that the van der Waals coefficients are very sensitive 
to the value of the dielectric constant. 



IV. MODELS 

Our proposed ab initio techniques are quite efficient 
and allow the calculation of van der Waals coefficients 
for systems with ~300 atoms, even in relatively mod- 
est computer systems. Moreover, the crystals of Si and 
SiC considered here contain only two atoms per unit cell 
and allow very fast calculations. Nevertheless, a full ab 
initio study of dispersion interactions of large nanostruc- 
tures/biological molecules on complex surfaces can be- 
come a computationally demanding task. Therefore, it is 
desirable to design accurate model van der Waals poten- 
tials, based on ab-initio calculations, to be used for these 
calculations, where a number of atoms of the order of 
1000 and even larger can be easily attained. Two prob- 
lems need to be addressed: i) how to model the dynamic 
polarizability at imaginary frequency of the nanoobjcct, 
ii) how to model the dielectric function at imaginary fre- 
quency of the solid. 



A. Model for the dielectric function 

Let us start by the modelling of the semiconductor or 
insulator surface. We have seen that the simplest expres- 
sion for the longitudinal dielectric function can be derived 
by applying Fermi's golden rule and assuming the case 
of a (nearly) homogeneous material: 



, 1 8tt I ^ |(' 
e q,w = 1+— 77 V — 

H Pi 1 



= iqr 



up 



irj 



(15) 

We have observed in Section ITlTI that the RPA approxi- 
mation gives a rather good dielectric function at imagi- 
nary frequency thanks to a compensation of errors. It is 
thus reasonable to start from this simple approximation 
to design an analytical model. A crude approximation 
consists in replacing the energy differences in Eq. (|15|) 
with some average excitation energy w av . By exploiting 



3 



12 
10 

8 
6 
4 



Si 


model 

;x a + (3 or/q 2 ■ ■ ■ 
NX exp. 


SiC 






0.1 1 




co (a.u.) 



FIG. 7: (Color online) Dielectric function of Si and SiC along 
the imaginary axis: the model function is compared with the 
curve obtained using the LRC dynamical approximation for 
the xc kernel, which gives the best agreement with the exper- 
iment. 



the sum rules and the asymptotic behavior of cm(w) one 
gets in the optical limit (q — > 0)^: 



emM = 1 - 



with 



H 



lrju 



(16) 



Equation (|16p has the same form as the Lorentz dielectric 
function of bound charged carriers with frequency uj av . 
This model predicts for the static dielectric constant 



cm(0) 



p 
,2 



(17) 



This equation can be used together with the value of the 
plasma frequency uj p = AirN c i/V to estimate cj av . In this 
approximation one only needs to fix two parameters that 
are easily accessible from experiments, the static dielec- 
tric constant and the volume of the unit cell. At imagi- 
nary frequencies an analogous expression can be derived 
as a function of the same parameters: 



e M (m) = 1 + 



(18) 



We remind that £m(iw) is a real function. The resulting 
model dielectric function is plotted in Fig. [7] for both Si 
and SiC, together with the best theoretical curve (the 
one obtained using the dynamical LRC model for the 
xc kernel) and the experimental curve for Si. The model 
function is substantially worse than the calculated curves, 
even with respect to the RPA calculations. 



B. Models for the cluster dynamical polarizabilities 

For the atomic polarizabilities at imaginary frequency 
a(iu) it is convenient to use to already mentioned London 
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TABLE I: Values for the static polarizability a(0), and for the Hamaker C6 coefficient. We show both the calculated values 
with (TD)DFT, as well as the values obtained using a bond polarization model (BPM) and effective medium theory (EMT), 
and the respective percent errors. 
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approximation, Eq. (|12[) . In that case, two parameters 
for every cluster have to be fixed: its static polarizabil- 
ity a(0) and the energy of an effective frequency u>\. Of 
course, one does not want to define a different set of pa- 
rameters for every cluster, but only two parameters for all 
possible clusters of a fixed species. In the case of a larger 
molecule or a cluster, the static polarizability can be es- 
timated by making use of the bond polarization model 
(BPM) 47 ' 48 as suggested by Jiemchooroj et al^ In this 
model, the total polarizability is obtained summing over 
the contributions from the individual polarizable entities: 
the covalent bonds. In our case there are only two kinds 
of bonds, Si-Si and Si-H, and we can write the static 
polarizability as 

Oi(0) = nf- Si a Si -si + ™f " H asi-H , (19) 

where nf 1_Sl and nf 1_H are the number of Si-Si and Si- 
ll bonds, respectively, of the cluster i. Here we have 
indicated with asi-Si and asm the contributions to the 
polarizability due to the Si-Si and the Si-H bonds, respec- 
tively. Upon substitution in the integral the London 
model gives for a homo-molecular Hamaker constant the 
expression (|13|) . and for the hetero- molecular Hamaker 



constant: 

Cff =3 ^ ai (0)a 3 -(0), (20) 

where (fl~9| gives the values of cti(0). If we perform a set of 
ab-initio calculations of Cg l and aj(0) for small- medium 
size clusters, we can extract u>\ from Eq. (fl~3"|) and the 
parameters asi-Sb c*Si-H fitting the theoretical curve 
for the static polarizability with Eq. (flUl) . 

The small dispersion of values for uj\ in Fig. [^suggests 
that it is possible to determine a single average frequency 
H>i =0.343 Ha for all clusters. For the small nanocrystals 
with less than 10 Si atoms this approximation is not very 
precise, but we are interested in getting information on 
the interaction of larger systems, that cannot be easily 
studied by ab-initio techniques. The parameters asi-Si 
and asm are fixed by fitting Eq. (|19[) to the curve for the 
static polarizability calculated within the TDLDA (see 
Fig. [2]). The outcome are the values asi-Si =13.41 a.u. 
and cvsi-H =10.56 a.u. In Table [I] we can verify the excel- 
lent agreement (|Aa(0)| < 1% for all the clusters, exclud- 
ing the smallest ones) between the static polarizabilities 
calculated using TDLDA and the additivity model for the 
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nanocrystals. The same agreement is conserved for the 
estimated values of Cq: For the big clusters, the differ- 
ence with the calculated values is remarkably small (see 
Table U. 

A different approach to model the dynamic polarizabil- 
ity of a nanocrystal is to start from the dielectric function 
of the corresponding bulk crystal and apply the effec- 
tive medium theory (EMT).— This classical approach is 
based on the solution of Maxwell's equations with the as- 
sumption that the dielectric response of each constituent 
of the system is the one of the corresponding bulk. This 
assumption is better justified when the size of the com- 
posing objects is large. In fact, EMT completely ne- 
glects the microscopic scale details, such as atoms and 
bonds. In this respect it is complementary to the additive 
procedure. However, it handles correctly the boundary 
conditions for the Maxwell's equations at the interfaces, 
which give very important contributions to the dielectric 
response through the crystal local field effects. Our clus- 
ters can be considered as a sphere of Si in vacuum with 
a filling factor / that goes to zero. The Maxwell-Garnett 
expression 51 yields in this specific case: 



3m {a (ui)} = 



- Jm 



4tt 



M + 2 



(21) 



where e|j is the complex dielectric function of bulk sili- 
con, and where V s is the volume of the spherical cluster. 
By applying the analogous of the Laplace transformation 
(TT4")) to the dynamical polarizability at real frequencies 
and using once again the single-oscillator model (| 16[) for 
the dielectric function of bulk Si, we obtain for the dy- 
namical polarizability at imaginary frequencies: 



a (iu) 



4tt u 2 + ojI v + cj2/3 ' 



(22) 



This expression can be rewritten in the same form as the 
London model by imposing: 



and 



a(0) 



V s 



47T wl 



■ W g/3 



^av + w p/3- 



(23) 



(24) 



With respect to the model for the dielectric function of 
the crystal we have here an extra parameter to be esti- 
mated: the volume of the spherical cluster. In the limit 
of a very large cluster, we can assume that the volume 
per Si atom is the same as in the bulk crystal. To ob- 
tain better results for small-medium sized clusters it is 
necessary to include the contribution to the volume due 
to the hydrogen atoms at the surface (considering a Si-H 
bond distance of about 1.5 A, we can assume a volume 
of about 11.4 a.u. per hydrogen atom). The value of u)\ 
is 0.4 Ha, which is 20% larger than u>\ evaluated in the 



BPM. The values of a(0), on the other hand, are system- 
atically smaller than their counterparts in the BPM. As 
a result, Cq coefficients are underestimated by about 5% 
for the larger clusters and up to 10-15% for the smaller 
ones. 

The non-trivial advantage of this second approach is 
the fact that no ab-initio calculation needs to be per- 
formed to fit a(0) and w\. 



C. Modelled Cj, coefficients 

Finally, we calculated the C3 coefficients both for Si 
nanocrystals on Si surfaces and on SiC surfaces com- 
bining the single-oscillator model for the surface the the 
bond polarization model for the nanocrystal. The curves 
obtained are plotted in Figs. [3] and 2] where they can be 
compared with the results of the ab-initio calculations. 
The model calculations are rigidly shifted to higher val- 
ues by about 5% for Si and 7% for SiC (when compared 
with the dynamical LRC model) . Curiously, they almost 
overlap with the results obtained using the LRC xc kernel 
for the determination of the dielectric function. This can 
be understood by inspecting Fig. [7] the error is likely to 
be entirely due to the insufficiently accurate description 
of the dielectric function of the bulk material. 



V. CONCLUSIONS 

We have demonstrated how the leading terms of the 
van der Waals forces acting between nanostructures, and 
between nanostructures and non-metallic surfaces, can be 
accurately and inexpensively computed from first prin- 
ciples. The key ingredients can be reliably obtained 
with state-of-the-art theoretical schemes and computa- 
tional procedures. The dynamical polarizabilities of large 
nanostructures at imaginary frequencies, for example, 
can be safely computed with the Sternheimer reformu- 
lation of time-dependent density-functional theory. Re- 
garding the other key ingredient necessary to obtain 
the cluster-surface interaction, the macroscopic dielec- 
tric constant, it can also be computed by making use 
of TDDFT. In this case, special care has to be taken 
with the choice of the xc kernel. We have found that a 
particularly simple and reliable scheme is to make use of 
the dynamical "long-range contribution" (LRC) kernel^ 
However, if the bulk to be studied lies outside the limits 
of validity of the dynamical LRC model (large gap insula- 
tors, systems with strongly bound excitons), one proba- 
bly has to resort to the full solution of the Bethe-Salpeter 
equation (or, equivalently, to a TDDFT calculation based 
on a kernel derived from the Bethe-Salpeter equation). 34 
This necessarily implies larger computational costs. 

We also suggest some simplified models that should 
supply reasonable estimates for the van der Waals coef- 
ficients, for those cases in which first principle calcula- 
tions are out of range. We have found that modelling 
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the bulk dielectric function can lead to a substantial er- 
ror (~ 10%). However, using a bond polarization model 
or the effective medium theory for the nanostructure dy- 
namical polarizability, yields very precise results, espe- 
cially for large systems. This opens the way to the simu- 
lation of the van der Waals interaction for large nanocrys- 
tals. 
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